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The iterated Crank-Nicolson is a predictor-corrector algorithm commonly used in numerical relativity for 
the solution of both hyperbolic and parabolic partial differential equations. We here extend the recent work on 
the stability of this scheme for hyperbolic equations by investigating the properties when the average between 
the predicted and corrected values is made with unequal weights and when the scheme is applied to a parabolic 
equation. We also propose a variant of the scheme in which the coefficients in the averages are swapped between 
two corrections leading to systematically larger amplification factors and to a smaller numerical dispersion. 
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I. INTRODUCTION 



resentation of the time derivative in eq. Q 



In a recent paper Ql, the stability of the iterated Crank- 
Nicolson (ICN) method was investigated for the solution of 
hyperbolic partial differential equations. We here extend the 
work in ref. 1 1] in three different ways. Firstly, we investigate 
the stability properties of the ICN method when the average 
between the predicted and corrected values is made with un- 
equal weights as recently used in refs. Secondly, we 
apply the above analysis to a prototypical parabolic partial dif- 
ferential equation, whose solution is also becoming important 
within numerical relativity simulations |3]. Finally, we pro- 
pose a variant of the scheme, valid for both hyperbolic and 
parabolic equations, in which the coefficients in the averages 
are swapped between two corrections, leading to larger ampli- 
fication factors and smaller numerical dispersion. 

The paper is organized as follows: in SectionsllJlandlllllwe 
recall the definition of the ICN as a predictor-corrector method 
and as a ^-method, respectively. In Sections HVl andlVl on the 
other hand, we discuss the stability properties of the 0-ICN in 
the case of hyperbolic and parabolic equations, respectively. 
The analysis of the truncation error, numerical dissipation and 
dispersion is presented in Section [V71 while the conclusions 
are collected in Section lVTll 

II. ICN AS A PREDICTOR-CORRECTOR METHOD 

Restricting our discussion to one spatial dimension, here- 
after we will consider a first-order in time partial differential 
equation of the type 



du{x, t) 
dt 



C{u{x, f)) 



(1) 



where £ is a generic quasi-linear partial differential operator 
which we assume to contain either first-order or second-order 
spatial partial derivatives. Most equations in numerical rela- 
tivity can be recast in this form and more complex operators 
follow from these two cases (see |4] and references therein). 

After introducing a discretization Ax in space and At in 
time, and truncating at the first order the finite-difference rep- 



du 
dt 



-T 1 ~ "J 

At 



+ O(At) , 



the generic solution of Q can be expressed as 



< +1 = »? 



AtL(ut) 



(2) 



(3) 



where, as usual, u" = u(jAx,nAt) with j and n integers, 
and L is the finite-difference form of the differential operator 
C. The spatial index k varies according to the order at which 
the operator is represented, with k = j ± 1 for a second- 
order accurate, first-order spatial derivative [cf. eq.(|5)], or 
with k = j, j ± 1 for a second-order accurate, second-order 
spatial derivative [cf. ea.( ll8H . 

The ICN scheme discussed in [ 1 ] is then the modification of 
the implicit Crank-Nicolson scheme [5] as obtained by trun- 
cating, at some point, the following infinite sequence of pre- 
dictions and corrections 



^u- +1 =u j+ AtL(^u^ 1/2 ) , 



^u] +1 = u] + AtL ( (2) < +1 



/2 



(4a) 
(4b) 
(4c) 
(4d) 
(4e) 



where Wu" +1/2 is the M-th average and ( M )«™ +1 , 
me ^f. m predicted and corrected solutions, re- 
spectively. 



III. ICN AS A 0-METHOD 

In the ICN method the M-th average is made weighting 
equally the newly predicted solution ( M )it™ +1 and the solution 



2 



at the "old" timelevel" u n . This, however, can be seen as the 
special case of a more generic averaging of the type 



(A/) -n+l/2 = Q (M)~n+1 + g y 



(5) 



where < < 1 is a constant coefficient. Predictor-corrector 
schemes using this type of averaging are part of a large class 
of algorithms named 9-methods |6], and we refer to the ICN 
generalized in this way as to the "0-ICN" method. 

A different and novel generalization of the 0-ICN can be 
obtained by swapping the averages between two subsequent 
corrector steps, so that in the M-th corrector step 

(M) -n+l/2 = (1 _ Q) (M)-„+l + 0u n ^ (fi) 

while in the (M + l)-th corrector step 

(Af+l) -n+l/2 = e (M+l)-„+l + (1 _ e)u n (?) 

Note that as long as the number of iterations is even, the se- 
quence in which the averages are computed is irrelevant. In- 
deed, the weights and 1 — 9 in eqs. could be inverted 
and all of the relations discussed hereafter for the swapped 
weighted averages would continue to hold after the transfor- 
mation 0^1 — 0. 

Although the properties of the 0-ICN do not seem to have 
been discussed before, the scheme has already found appli- 
cation in numerical relativity calculations, where it has been 
used with a coefficient 9 = 0.6 in the solution of the rela- 
tivistic hydrodynamics equations for ideal |2] and viscous flu- 
ids 1 3]. In these works it was found that the use of a weighting 
coefficient different from 1/2 yielded "an improved stability". 
In Sect. I VII we will show that such a choice has effectively 
only increased the numerical dissipation of the scheme. 

IV. HYPERBOLIC EQUATIONS 

To discuss the properties of the 0-ICN we consider as model 
hyperbolic equation the one-dimensional advection equation 



du du 
dt dx 



(8) 



where v is a constant coefficient. A second-order accu- 
rate finite-difference representation of the right-hand-side of 
eq. (|8} is then easy to derive and has the form 



2Ax 



0(Ax 



2^ 



(9) 



A. Constant Arithmetic Averages 

Using a von Neumann stability analysis, Teukolsky has 
shown that for a hyperbolic equation the ICN scheme with 
M iterations has an amplification factor 1 1 ] 



M 



(M) 



(10) 



where (3 = v[At/(2Ax)] sin kAx d. We recall that in a 
von Neumann stability analysis the eigenmodes of the finite- 
difference equations are expressed as u™ = £ n e lk i Ax ^ where 
k is a real spatial wavenumber and £ = £(fc) is a complex 
number. Stability then requires that |£| = \/t[£* < 1 and 
in eq. ( 1101 this leads to an alternating pattern in the number 
of iterations. More specifically, zero and one iterations yield 
an unconditionally unstable scheme, while two and three it- 
erations a stable one provided that (3 2 < 1; four and five 
iterations lead again to an unstable scheme and so on. Fur- 
thermore, because the scheme is second-order accurate from 
the first iteration on, Teukolsky's suggestion when using the 
ICN method for hyperbolic equations was that two iterations 
should be used and no more d] . This is the number of itera- 
tions we will consider hereafter. 



B. Constant Weighted Averages 

Performing the same stability analysis for a 0-ICN is only 
slightly more complicated and truncating at two iterations the 
amplification factor is found to be 



£ = 1 - 2i/3 - 4(3 2 9 + 8i/3 



3/)2 



(ID 



where £ is a shorthand for ( 2 )£. The stability condition in this 
case translates into requiring that 



16/? 4 4 - 4/? 2 2 - 29 + 1 < 
or, equivalently, that for 9 > 3/8 



(12) 



20 



20 



< (3 < 



29 



29 



(13) 



which reduces to (P < 1 when = 1/2. Because the con- 
dition ( II 31 must hold for every wavenumber k, we consider 
hereafter j3 = vAt/ (2 Ax) and show in the left panel of Fig. ^ 
the region of stability in the (0, (3) plane. The thick solid lines 
mark the limit at which |£| = 1, while the dotted contours 
indicate the different values of the amplification factor in the 
stable region. 

A number of comments are worth making. Firstly, although 
the condition ( 1131 allows for weighting coefficients < 1/2, 
the 0-ICN is stable only if > 1/2. This is a known property 
of the weighted Crank-Nicolson scheme |6] and inherited by 
the 0-ICN. In essence, when 0^1/2 spurious solutions ap- 
pear in the method [7] and these solutions are linearly unstable 
if < 1/2, while they are stable for > 1/2 01 (An alter- 
native and simpler explanation is also presented in Sect. IVIi . 
For this reason we have shaded the area with < 1/2 in the 
left panel of Fig.^to exclude it from the stability region. Sec- 
ondly, the use of a weighting coefficient > 1/2 will still lead 
to a stable scheme provided that the timestep (i.e., (3) is suit- 
ably decreased. Finally, as the contour lines in the left panel 
of Fig. \l\ clearly show, the amplification factor can be very 
sensitive on 0. 
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FIG. 1: Left panel: stability region in the (0,/J) plane for the two-iterations 6-lCN for the advection equation 0. Thick solid lines mark 
the limit at which |£| = 1, while the dotted contours indicate the values of the amplification factor in the stable region. The shaded area for 
9 < 1/2 refers to solutions that are linearly unstable 1 8]. Right panel: same as in the left panel but when the averages between two corrections 
are swapped. Note that the amplification factor in this case is less sensitive on 9 and always larger than the corresponding amplification factor 
in the left panel. 



C. Swapped weighted averages 



V. PARABOLIC EQUATIONS 



The calculation of the stability of the (9-ICN when the 
weighted averages are swapped as in eqs. l|6} and @ is some- 
what more involved; after some lengthy but straightforward 
algebra we find the amplification factor to be 



£ = 1 - 2i/3 - 4/3 2 (9 + 8i/3 3 (9(l 



(14) 



which differs from il It only in that the 9 2 coefficient of the 
0(P 3 ) term is replaced by 8(1 — 8). The stability requirement 
|£| < 1 is now expressed as 



16/3 4 2 (1 



- 4/3 2 6»(2 - 38) - 28 + 1 < 



(15) 



Solving the condition ( fTSl with respect to (5 amounts then to 
requiring that 



(3> 



0< 



V2-38-V48-118 2 



2(1 



2fV 



\j2 - 38 + V48 - 118' 2 



2(1 -8)V28 



(16a) 
(16b) 



which is again equivalent to f3 2 < 1 when 8 — 1/2. The 
corresponding region of stability is shown in right panel of 
Fig. ^ an d should be compared with left panel of the same 
Figure. Note that the average-swapping has now considerably 
increased the amplification factor, which is always larger than 
the corresponding one for the 6>-ICN in the relevant region of 
stability (i.e., for 1/2 < 6 < 1 fil l. 



We next extend the stability analysis of the 0-ICN to the a 
parabolic partial differential equation and use as model equa- 
tion the one-dimensional diffusion equation 



du £)® 2u Q 
dt dx 2 



(17) 



where D is a constant coefficient which must be positive for 
the equation to be well-posed. 

Parabolic equations are commonly solved using implicit 
methods such as the Crank-Nicolson, which is uncondition- 
ally stable and thus removes the constraints on the timestep 
[i.e., At ps 0(Ax 2 )] imposed by explicit schemes j^]. 
In multidimensional calculations, however, or when the set 
of equations is of mixed hyperbolic-parabolic type, implicit 
schemes can be cumbersome to implement since the resulting 
system of algebraic equations does no longer have simple and 
tridiagonal matrices of coefficients. In this case, the most con- 
veniente choice may be to use an explicit method such as the 
ICN. 

Also in this case, the first step in our analysis is the deriva- 
tion of a finite-difference representation of the right-hand-side 
of eq. illl which, at second-order, has the form 



2u° 



Ax 2 



0(Ax 2 ) 



(18) 
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A. Constant Arithmetic Averages 



C. Swapped Weighted Averages 



Next, we consider first the case with constant arithmetic av- 
erages (i.e., 9 = 1/2) and the expression for the amplification 
factor after M-iterations is then purely real and given by 



M 



(M) 



£=l + 2£(-7r 



(19) 



where 7 = (2DAt/Ax 2 ) sin 2 (fcAa;/2). Requiring now for 
stability that yj^ 2 < 1 ar, d bearing in mind that 



M 



-l<^(-7)" +1 <0, for 7 <l, 



(20) 



n=0 



we find that the scheme is stable for any number of iterations 
provided that 7 < 1. Furthermore, because the scheme is 
second-order accurate from the first iteration on, our sugges- 
tion when using the ICN method for parabolic equations is 
that one iteration should be used and no more. In this case, in 
particular, the ICN method coincides with a FTCS scheme 1 9 ] . 

Note that the stability condition 7 < 1 introduces again a 
constraint on the timestep that must be At < Air 2 /(2D) and 
thus 0(Ax 2 ). As a result and at least in this respect, the ICN 
method does not seem to offer any advantage over other ex- 
plicit methods for the solution of a parabolic equation [14]. 



After some lengthy algebra the calculation of the amplifi- 
cation factor for the 0-ICN method with swapped weighted 
averages yields 



£ = 1 - 27 + 4 7 2 6> - 8 7 3 0(1 - 



and stability is then given by 

-1 < 1- 2 7 + 4 7 2 6 



87 3 6>(1 - 9) < 1 



(24) 



(25) 



Note that none of the two inequalities is always true and in 
order to obtain analytical expressions for the stable region we 
solve the condition d25l with respect to 9 and obtain 



< 



< 



> 



27 - 1 + y/^f -47 + 5 
4^ ' 

7(27 - 1) - ^7 (47 3 - 47 2 + 57- 4) 
4^2 



7(27 - 1) + 1/7 (4 7 3 - 47 2 + 57-4) 
I7 1 



(26a) 
(26b) 
(26c) 



The resulting stable region for sin kAx = 1 is plotted in the 
right panel of Fig.|2]and seems to suggest that arbitrarily large 
values of 7 could be considered when 9 > 0.6 It should be 
noted, however, that the amplification factor is also severely 
reduced as larger values of 7 are used and indeed it is essen- 
tially zero in the limit 9 —> 1. 



B. Constant Weighted Averages 



VI. TRUNCATION ERROR, DISSIPATION AND 
DISPERSION 



We next consider the stability of the 0-ICN method but fo- 
cus our attention on a two-iterations scheme since this is the 
number of iterations needed in the solution of the parabolic 
part in a mixed hyperbolic -parabolic equation when, for in- 
stance, operator-splitting techniques are adopted [9]. In this 
case, the amplification factor is again purely real and given by 



3/)2 



£ = 1 - 2 7 + 4 7 ^6» - 87 J 6> 
so that stability is achieved if 

< 7 (1 - 2(9 7 + 46>V) < 1 



(21) 



(22) 



Since 7 > by definition, the left inequality is always satis- 
fied, while the right one is hue provided that, for 7 < 4/3, 



7- V7(4~3 7 ) 
4 7 2 



< 9 < 



7+ y/ 7 (4-3 7 ) 
4 7 2 



(23) 



The stability region described by the condition (I23> is 
shown in the left panel of Fig.|2]for sinfcAx = 1 and illus- 
trates that the scheme is stable for any value < 9 < 1, and 
also that slightly larger timesteps can be taken when 9 ~ 0.2. 



Although not often appreciated, the 0-ICN method is only 
first-order accurate in time as an obvious consequence of the 
first-order approximation in the time derivative [cf. eq. 0]. 
However, this is not hue if 9 = 1/2, in which case the method 
becomes second-order in both space and time. 

To appreciate this in the case of the advection equation 0, 
we report the finite-difference expressions for the time and 
spatial derivatives in eq. (|8j, writing out explicitly the coeffi- 
cients of the O(At) and 0(Ax 2 ) terms 



At 



du 



2Aa 



du 
d.r 



l a u 
2W 



At + 0(At z ) , (27) 



ld 3 u 
6 dx 3 



Ax 2 +0{Ax 4 ) 

n,j 



(28) 



The resulting local truncation error is then 



1 

2 ~ 

„.3fl2 



d 2 u 



dxdt 



At 



v d 3 u 
6 dx 3 



Ax 2 



d 3 u 
dx 3 



At 2 - 



O (At 3 , A 



"■:) 

3'\ 



1 d 3 u 

6 at 3 " 



At 2 



(29) 
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FIG. 2: Left panel: stability region in the (8, 7) plane for the two-iterations 0-ICN for the diffusion equation i 171 . Thick solid lines mark the 
limit at which £ 2 = 1, while the dotted contours indicate the values of the amplification factor in the stable region. Right panel: same as in the 
left panel but with swapping the averages between two corrections. 



clearly indicating that the (9-ICN is generally only first-order 
accurate in time, becoming second-order if 6 = 1/2. The 
truncation error is also useful to quantify the numerical dis- 
sipation and dispersion inherent to the 0-ICN method. Using 
eq. (|8} to replace the time derivative with a spatial one, in 
fact, eq. ( I29> shows that the 0-ICN introduces a dissipative 
term proportional to d 2 u/dx 2 and with coefficient 



(30) 



In other words, the 0-ICN is intrinsically dissipative, with 
a dissipation coefficient that is generically O(At) and 
0(At 3 , Ax 3 ) only when 9 = 1/2. Furthermore, it is now 
apparent why 9 must be larger or equal to 1/2; any choice dif- 
ferent from this, in fact, would change the sign of e a d v , leading 
to an ill-posed equation with exponentially growing solutions 
[cf. eq. Ol- 

Expression (I30i also clarifies the behaviour found in 
refs. ||2j, |3J. Since stability in a numerical scheme is either 
gained or lost but cannot be "improved", the use of a weight- 
ing coefficient 9 > 1/2 (and of a suitable timestep) has sim- 
ply the effect of increasing the numerical dissipation of the 
scheme. Of course, this is often a desirable feature to suppress 
the growth of instabilities, as in the case of the Lax-Friedrichs 
scheme, whose numerical dissipation stabilizes the otherwise 
unconditionally unstable FTCS scheme |9]. 

An alternative route to a second-order, moderately dissipa- 
tive scheme is to choose 



Ax 2 
vAt 



(31) 



with Ax 2 / (vAt) < 1/2, so that the leading error-term in i29i 
becomes again 0(At 2 , Ax 2 ). A prescription of the type ( 13 U 
may be the optimal one for the 9-lCN method as it provides a 
small amount of numerical dissipation and reduces the trun- 
cation error. 

The truncation error i29\ also indicates that the 0-ICN in- 
troduces a dispersive term proportional to d 3 u/dx 3 given by 



1 

Xadv ~~ \ g 



1 



24/3 2 



v 3 At 2 



(32) 



and responsible, for instance, for different propagation speeds 
of the Fourier modes in the initial data {i.e., phase drifts). 

All what discussed so far for the 6>-ICN scheme continues to 
hold also when the averages are swapped, the only difference 
being that the dispersive contribution is instead given by 



Xadv 



1 



24/3 5 



■/At 2 



(33) 



and is therefore smaller for 9 > 1/2, making this variant to 
the 9-lCN preferable overall. 

Similar calculations can be carried out also for the parabolic 
equation dl7> and the local truncation error in this case is 



2 J dtdx 2 



At + (At 2 , Ax 2 



(34) 



indicating that mathematically the 0-ICN is again only first- 
order accurate in time, with second-order accuracy being 
recovered for 9 = 1/2. However, stability requires that 
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FIG. 3: Upper panel: Solution of the advection equation J5J us- 
ing the #-ICN method and shown after 10 crossing times. Different 
curves refer to either the analytic solution or to the numerical ones 
with different weighting coefficients. The small inset, instead, shows 
the smaller diffusion and dispersion obtained when the averages are 
swapped (see main text for details). Lower panel: L2 norms of the 
solutions in the upper panel plotted as a function of time. 

At = O (Arc 2 ) (cf. Sect. IV Al l and thus the truncation error 
is effectively e T = O (Ax 2 ) for all of the allowed values 
< B < 1. 

Finally, using eq. dl7> to the replace the time derivative 
in (I34> shows that the 0-ICN for a parabolic equation has an 
additional dissipative term proportional to d 4 u/dx 4 with co- 
efficient 

eaiff = r - 5) D2At > (35) 

which is again zero only for 9 = 1/2. 

As a purely representative example we show in Fig. [5] the 
application of the 0-ICN method for the solution of the ad- 
vection equation with v = 1 and (3 = 0.6 (cf. Fig. [Q. 
The numerical domain has length 1.0 and was covered with 
200 equally spaced gridpoints. The initial solution, given by 
a Gaussian centred at x = 0.5 and with variance 0.1, was 
evolved for 10 crossing times using periodic boundary condi- 
tions. Different curves in the upper panel refer to either the 
analytic solution at the final time (dotted line) or to the nu- 
merical solutions as obtained with different weighting coeffi- 
cients. Note that already with 9 = 1/2 (solid line) the numer- 
ical solution is slightly diffused but suffers from considerable 
dispersion as apparent from the considerable "delay" and the 



presence of negative values to the left of the maximum. These 
dispersion errors can be reduced if larger values of the weight- 
ing coefficients are used as indicated by the short-dashed and 
long-dashed lines referring to 9 = 0.6 and 9 = 0.8, respec- 
tively. This improvement, however, also comes with a larger 
dissipation and truncation error (as mentioned in Sect. lVIl the 
system is just first-order in time with 9 > 1/2) fl5ll . This is 
particularly evident when considering the evolution of the L2 
norms of the solutions as reported in the lower panel of Fig. [3] 
It is interesting to note that for 9 = 0.8 the L2 norm of the so- 
lution has been reduced of about 25% after 10 crossing times, 
while this decrease is less than 1% when 9 = 1/2. 

Finally, the two small insets in Fig.|3]offer a comparison in 
the solutions for 9 = 0.6 when the coefficients in the averages 
are either held constant (short-dashed lines) or swapped be- 
tween two subsequent corrector steps (dot-dashed lines). Al- 
though the difference is rather small for the selected set of pa- 
rameters, it is evident that the swapping of the coefficients has 
the effect of decreasing both the dispersion (the dot-dashed 
line in the upper inset has a smaller "delay") and the diffusion 
(at any given time the dot-dashed line in the lower inset has a 
larger value). 



VII. CONCLUSIONS 

We have extended the recent work on the properties of the 
ICN scheme for hyperbolic equations by investigating the sta- 
bility properties when it is treated as a 0-method, i.e., when the 
average between the predicted and corrected values is made 
with unequal weights. In addition we have studied the proper- 
ties of the 0-ICN method for a model parabolic equation and 
proposed a variant of the scheme, valid for both hyperbolic 
and parabolic equations, in which the unequal coefficients co- 
efficients in the averages are swapped between two subsequent 
corrector steps. This novel approach leads to amplification 
factors that are systematically larger than those found in the 
6 -ICN method and to a smaller numerical dispersion. 

Overall, our results indicate that although generally only 
first-order accurate in time, the 0-ICN method is a flexible ap- 
proach to the time-integration of partial differential equations, 
particularly when these are of mixed hyperbolic-parabolic 
type. Because the use of unequal coefficients in the average 
provides a small but nonzero amount of numerical dissipa- 
tion, this could prove useful in numerical relativity calcula- 
tions which may suffer from the development of numerical 
instabilities and for which lower-order evolution schemes are 
an acceptable compromise between accuracy and stability. 
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